home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / AA_51.ZIP / NUTATE.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  10KB  |  393 lines
  1.  
  2. /* Nutation in longitude and obliquity
  3.  * computed at Julian date J.
  4.  *
  5.  * References:
  6.  * "Summary of 1980 IAU Theory of Nutation (Final Report of the
  7.  * IAU Working Group on Nutation)", P. K. Seidelmann et al., in
  8.  * Transactions of the IAU Vol. XVIII A, Reports on Astronomy,
  9.  * P. A. Wayman, ed.; D. Reidel Pub. Co., 1982.
  10.  *
  11.  * "Nutation and the Earth's Rotation",
  12.  * I.A.U. Symposium No. 78, May, 1977, page 256.
  13.  * I.A.U., 1980.
  14.  *
  15.  * Woolard, E.W., "A redevelopment of the theory of nutation",
  16.  * The Astronomical Journal, 58, 1-3 (1953).
  17.  *
  18.  * This program implements all of the 1980 IAU nutation series.
  19.  * Results checked at 100 points against the 1986 AA; all agreed.
  20.  *
  21.  *
  22.  * - S. L. Moshier, November 1987
  23.  * October, 1992 - typo fixed in nutation matrix
  24.  */
  25.  
  26. /* The answers are posted here by nutlo():
  27.  */
  28. double jdnut = -1.0;    /* time to which the nutation applies */
  29. double nutl = 0.0;    /* nutation in longitude (radians) */
  30. double nuto = 0.0;    /* nutation in obliquity (radians) */
  31. extern double eps, coseps, sineps, DTR;
  32.  
  33. /* Each term in the expansion has a trigonometric
  34.  * argument given by
  35.  *   W = i*MM + j*MS + k*FF + l*DD + m*OM
  36.  * where the variables are defined below.
  37.  * The nutation in longitude is a sum of terms of the
  38.  * form (a + bT) * sin(W). The terms for nutation in obliquity
  39.  * are of the form (c + dT) * cos(W).  The coefficients
  40.  * are arranged in the tabulation as follows:
  41.  * 
  42.  * Coefficient:
  43.  * i  j  k  l  m      a      b      c     d
  44.  * 0, 0, 0, 0, 1, -171996, -1742, 92025, 89,
  45.  * The first line of the table, above, is done separately
  46.  * since two of the values do not fit into 16 bit integers.
  47.  * The values a and c are arc seconds times 10000.  b and d
  48.  * are arc seconds per Julian century times 100000.  i through m
  49.  * are integers.  See the program for interpretation of MM, MS,
  50.  * etc., which are mean orbital elements of the Sun and Moon.
  51.  *
  52.  * If terms with coefficient less than X are omitted, the peak
  53.  * errors will be:
  54.  *
  55.  *   omit    error,          omit    error,
  56.  *   a <    longitude      c <    obliquity
  57.  * .0005"    .0100"        .0008"    .0094"
  58.  * .0046    .0492        .0095    .0481
  59.  * .0123    .0880        .0224    .0905
  60.  * .0386    .1808        .0895    .1129
  61.  */
  62. short nt[105*9] = {
  63.  0, 0, 0, 0, 2, 2062, 2,-895, 5,
  64. -2, 0, 2, 0, 1, 46, 0,-24, 0,
  65.  2, 0,-2, 0, 0, 11, 0, 0, 0,
  66. -2, 0, 2, 0, 2,-3, 0, 1, 0,
  67.  1,-1, 0,-1, 0,-3, 0, 0, 0,
  68.  0,-2, 2,-2, 1,-2, 0, 1, 0,
  69.  2, 0,-2, 0, 1, 1, 0, 0, 0,
  70.  0, 0, 2,-2, 2,-13187,-16, 5736,-31,
  71.  0, 1, 0, 0, 0, 1426,-34, 54,-1,
  72.  0, 1, 2,-2, 2,-517, 12, 224,-6,
  73.  0,-1, 2,-2, 2, 217,-5,-95, 3,
  74.  0, 0, 2,-2, 1, 129, 1,-70, 0,
  75.  2, 0, 0,-2, 0, 48, 0, 1, 0,
  76.  0, 0, 2,-2, 0,-22, 0, 0, 0,
  77.  0, 2, 0, 0, 0, 17,-1, 0, 0,
  78.  0, 1, 0, 0, 1,-15, 0, 9, 0,
  79.  0, 2, 2,-2, 2,-16, 1, 7, 0,
  80.  0,-1, 0, 0, 1,-12, 0, 6, 0,
  81. -2, 0, 0, 2, 1,-6, 0, 3, 0,
  82.  0,-1, 2,-2, 1,-5, 0, 3, 0,
  83.  2, 0, 0,-2, 1, 4, 0,-2, 0,
  84.  0, 1, 2,-2, 1, 4, 0,-2, 0,
  85.  1, 0, 0,-1, 0,-4, 0, 0, 0,
  86.  2, 1, 0,-2, 0, 1, 0, 0, 0,
  87.  0, 0,-2, 2, 1, 1, 0, 0, 0,
  88.  0, 1,-2, 2, 0,-1, 0, 0, 0,
  89.  0, 1, 0, 0, 2, 1, 0, 0, 0,
  90. -1, 0, 0, 1, 1, 1, 0, 0, 0,
  91.  0, 1, 2,-2, 0,-1, 0, 0, 0,
  92.  0, 0, 2, 0, 2,-2274,-2, 977,-5,
  93.  1, 0, 0, 0, 0, 712, 1,-7, 0,
  94.  0, 0, 2, 0, 1,-386,-4, 200, 0,
  95.  1, 0, 2, 0, 2,-301, 0, 129,-1,
  96.  1, 0, 0,-2, 0,-158, 0,-1, 0,
  97. -1, 0, 2, 0, 2, 123, 0,-53, 0,
  98.  0, 0, 0, 2, 0, 63, 0,-2, 0,
  99.  1, 0, 0, 0, 1, 63, 1,-33, 0,
  100. -1, 0, 0, 0, 1,-58,-1, 32, 0,
  101. -1, 0, 2, 2, 2,-59, 0, 26, 0,
  102.  1, 0, 2, 0, 1,-51, 0, 27, 0,
  103.  0, 0, 2, 2, 2,-38, 0, 16, 0,
  104.  2, 0, 0, 0, 0, 29, 0,-1, 0,
  105.  1, 0, 2,-2, 2, 29, 0,-12, 0,
  106.  2, 0, 2, 0, 2,-31, 0, 13, 0,
  107.  0, 0, 2, 0, 0, 26, 0,-1, 0,
  108. -1, 0, 2, 0, 1, 21, 0,-10, 0,
  109. -1, 0, 0, 2, 1, 16, 0,-8, 0,
  110.  1, 0, 0,-2, 1,-13, 0, 7, 0,
  111. -1, 0, 2, 2, 1,-10, 0, 5, 0,
  112.  1, 1, 0,-2, 0,-7, 0, 0, 0,
  113.  0, 1, 2, 0, 2, 7, 0,-3, 0,
  114.  0,-1, 2, 0, 2,-7, 0, 3, 0,
  115.  1, 0, 2, 2, 2,-8, 0, 3, 0,
  116.  1, 0, 0, 2, 0, 6, 0, 0, 0,
  117.  2, 0, 2,-2, 2, 6, 0,-3, 0,
  118.  0, 0, 0, 2, 1,-6, 0, 3, 0,
  119.  0, 0, 2, 2, 1,-7, 0, 3, 0,
  120.  1, 0, 2,-2, 1, 6, 0,-3, 0,
  121.  0, 0, 0,-2, 1,-5, 0, 3, 0,
  122.  1,-1, 0, 0, 0, 5, 0, 0, 0,
  123.  2, 0, 2, 0, 1,-5, 0, 3, 0, 
  124.  0, 1, 0,-2, 0,-4, 0, 0, 0,
  125.  1, 0,-2, 0, 0, 4, 0, 0, 0,
  126.  0, 0, 0, 1, 0,-4, 0, 0, 0,
  127.  1, 1, 0, 0, 0,-3, 0, 0, 0,
  128.  1, 0, 2, 0, 0, 3, 0, 0, 0,
  129.  1,-1, 2, 0, 2,-3, 0, 1, 0,
  130. -1,-1, 2, 2, 2,-3, 0, 1, 0,
  131. -2, 0, 0, 0, 1,-2, 0, 1, 0,
  132.  3, 0, 2, 0, 2,-3, 0, 1, 0,
  133.  0,-1, 2, 2, 2,-3, 0, 1, 0,
  134.  1, 1, 2, 0, 2, 2, 0,-1, 0,
  135. -1, 0, 2,-2, 1,-2, 0, 1, 0,
  136.  2, 0, 0, 0, 1, 2, 0,-1, 0,
  137.  1, 0, 0, 0, 2,-2, 0, 1, 0,
  138.  3, 0, 0, 0, 0, 2, 0, 0, 0,
  139.  0, 0, 2, 1, 2, 2, 0,-1, 0,
  140. -1, 0, 0, 0, 2, 1, 0,-1, 0,
  141.  1, 0, 0,-4, 0,-1, 0, 0, 0,
  142. -2, 0, 2, 2, 2, 1, 0,-1, 0,
  143. -1, 0, 2, 4, 2,-2, 0, 1, 0,
  144.  2, 0, 0,-4, 0,-1, 0, 0, 0,
  145.  1, 1, 2,-2, 2, 1, 0,-1, 0,
  146.  1, 0, 2, 2, 1,-1, 0, 1, 0,
  147. -2, 0, 2, 4, 2,-1, 0, 1, 0,
  148. -1, 0, 4, 0, 2, 1, 0, 0, 0,
  149.  1,-1, 0,-2, 0, 1, 0, 0, 0,
  150.  2, 0, 2,-2, 1, 1, 0,-1, 0,
  151.  2, 0, 2, 2, 2,-1, 0, 0, 0,
  152.  1, 0, 0, 2, 1,-1, 0, 0, 0,
  153.  0, 0, 4,-2, 2, 1, 0, 0, 0,
  154.  3, 0, 2,-2, 2, 1, 0, 0, 0,
  155.  1, 0, 2,-2, 0,-1, 0, 0, 0,
  156.  0, 1, 2, 0, 1, 1, 0, 0, 0,
  157. -1,-1, 0, 2, 1, 1, 0, 0, 0,
  158.  0, 0,-2, 0, 1,-1, 0, 0, 0,
  159.  0, 0, 2,-1, 2,-1, 0, 0, 0,
  160.  0, 1, 0, 2, 0,-1, 0, 0, 0,
  161.  1, 0,-2,-2, 0,-1, 0, 0, 0,
  162.  0,-1, 2, 0, 1,-1, 0, 0, 0,
  163.  1, 1, 0,-2, 1,-1, 0, 0, 0,
  164.  1, 0,-2, 2, 0,-1, 0, 0, 0,
  165.  2, 0, 0, 2, 0, 1, 0, 0, 0,
  166.  0, 0, 2, 4, 2,-1, 0, 0, 0,
  167.  0, 1, 0, 1, 0, 1, 0, 0, 0,
  168. };
  169.  
  170.  
  171. /* arrays to hold sines and cosines of multiple angles
  172.  */
  173.  
  174. double ss[5][8];
  175. double cc[5][8];
  176. extern double ss[5][8], cc[5][8];
  177. int sscc(), epsiln(), showcor();
  178.  
  179. int nutlo(J)
  180. double J;
  181. {
  182. double f, g, T;
  183. double MM, MS, FF, DD, OM;
  184. double cu, su, cv, sv, sw;
  185. double C, D;
  186. int i, j, k, k1, m;
  187. short *p;
  188. double sin(), cos(), mod360();
  189.  
  190. if( jdnut == J )
  191.     return(0);
  192. jdnut = J;
  193.  
  194. /* Julian centuries from 2000 January 1.5,
  195.  * barycentric dynamical time
  196.  */
  197. T = (J-2451545.0)/36525.0;
  198.  
  199. /* Fundamental arguments in the FK5 reference system.
  200.  * The coefficients, originally given to 0.001",
  201.  * are converted here to degrees.
  202.  */
  203.  
  204. /* longitude of the mean ascending node of the lunar orbit
  205.  * on the ecliptic, measured from the mean equinox of date
  206.  */
  207. OM = ((2.22e-6*T + 0.00207833)*T - 1934.1362608)*T + 125.0445222;
  208. OM = DTR * mod360(OM);
  209.  
  210. /* mean longitude of the Sun minus the
  211.  * mean longitude of the Sun's perigee
  212.  */
  213. MS = ((-3.33e-6*T - 0.0001603)*T + 35999.0503400)*T + 357.5277233;
  214. MS = DTR * mod360(MS);
  215.  
  216. /* mean longitude of the Moon minus the
  217.  * mean longitude of the Moon's perigee
  218.  */
  219. MM = ((1.778e-5*T + 0.0086972)*T + 477198.8673981)*T + 134.9629814;
  220. MM = DTR * mod360(MM);
  221.  
  222. /* mean longitude of the Moon minus the
  223.  * mean longitude of the Moon's node
  224.  */
  225. FF  = ((3.056e-6*T - 0.00368250)*T + 483202.0175381)*T + 93.2719103;
  226. FF = DTR * mod360(FF);
  227.  
  228. /* mean elongation of the Moon from the Sun.
  229.  */
  230. DD  = (( 5.278e-6*T - 0.0019142)*T + 445267.1114800)*T + 297.8503631;
  231. DD = DTR * mod360(DD);
  232.  
  233. /* Calculate sin( i*MM ), etc. for needed multiple angles
  234.  */
  235. sscc( 0, MM, 3 );
  236. sscc( 1, MS, 2 );
  237. sscc( 2, FF, 4 );
  238. sscc( 3, DD, 4 );
  239. sscc( 4, OM, 2 );
  240.  
  241. /* first terms, not in table: */
  242. C = (-0.01742*T - 17.1996)*ss[4][0];    /* sin(OM) */
  243. D = ( 0.00089*T +  9.2025)*cc[4][0];    /* cos(OM) */
  244.  
  245. p = &nt[0]; /* point to start of table */
  246.  
  247. for( i=0; i<105; i++ )
  248.     {
  249. /* argument of sine and cosine */
  250.     k1 = 0;
  251.     cv = 0.0;
  252.     sv = 0.0;
  253.     for( m=0; m<5; m++ )
  254.         {
  255.         j = *p++;
  256.         if( j )
  257.             {
  258.             k = j;
  259.             if( j < 0 )
  260.                 k = -k;
  261.             su = ss[m][k-1]; /* sin(k*angle) */
  262.             if( j < 0 )
  263.                 su = -su;
  264.             cu = cc[m][k-1];
  265.             if( k1 == 0 )
  266.                 { /* set first angle */
  267.                 sv = su;
  268.                 cv = cu;
  269.                 k1 = 1;
  270.                 }
  271.             else
  272.                 { /* combine angles */
  273.                 sw = su*cv + cu*sv;
  274.                 cv = cu*cv - su*sv;
  275.                 sv = sw;
  276.                 }
  277.             }
  278.         }
  279. /* longitude coefficient */
  280.     f  = *p++ * 0.0001;
  281.     if( (k = *p++) != 0 )
  282.         f += 0.00001 * T * k;
  283.  
  284. /* obliquity coefficient */
  285.     g = *p++ * 0.0001;
  286.     if( (k = *p++) != 0 )
  287.         g += 0.00001 * T * k;
  288.  
  289. /* accumulate the terms */
  290.     C += f * sv;
  291.     D += g * cv;
  292.     }
  293. /*
  294. printf( "nutation: in longitude %.3f\", in obliquity %.3f\"\n", C, D );
  295. */
  296. /* Save answers, expressed in radians */
  297. nutl = DTR * C / 3600.0;
  298. nuto = DTR * D / 3600.0;
  299. return(0);
  300. }
  301.  
  302.  
  303.  
  304. /* Nutation -- AA page B20
  305.  * using nutation in longitude and obliquity from nutlo()
  306.  * and obliquity of the ecliptic from epsiln()
  307.  * both calculated for Julian date J.
  308.  *
  309.  * p[] = equatorial rectangular position vector of object for
  310.  * mean ecliptic and equinox of date.
  311.  */
  312. int nutate( J, p )
  313. double J;
  314. double p[];
  315. {
  316. double ce, se, cl, sl, sino, f;
  317. double dp[3], p1[3];
  318. int i;
  319. double sin(), cos();
  320.  
  321. nutlo(J); /* be sure we calculated nutl and nuto */
  322. epsiln(J); /* and also the obliquity of date */
  323.  
  324. f = eps + nuto;
  325. ce = cos( f );
  326. se = sin( f );
  327. sino = sin(nuto);
  328. cl = cos( nutl );
  329. sl = sin( nutl );
  330.  
  331. /* Apply adjustment
  332.  * to equatorial rectangular coordinates of object.
  333.  *
  334.  * This is a composite of three rotations: rotate about x axis
  335.  * to ecliptic of date; rotate about new z axis by the nutation
  336.  * in longitude; rotate about new x axis back to equator of date
  337.  * plus nutation in obliquity.
  338.  */
  339. p1[0] =   cl*p[0]
  340.         - sl*coseps*p[1]
  341.         - sl*sineps*p[2];
  342.  
  343. p1[1] =   sl*ce*p[0]
  344.         + ( cl*coseps*ce + sineps*se )*p[1]
  345.         - ( sino + (1.0-cl)*sineps*ce )*p[2];
  346.  
  347. p1[2] =   sl*se*p[0]
  348.         + ( sino + (cl-1.0)*se*coseps )*p[1]
  349.         + ( cl*sineps*se + coseps*ce )*p[2];
  350.  
  351. for( i=0; i<3; i++ )
  352.     dp[i] = p1[i] - p[i];
  353. showcor( "nutation", p, dp );
  354.  
  355. for( i=0; i<3; i++ )
  356.     p[i] = p1[i];
  357. return(0);
  358. }
  359.  
  360.  
  361. /* Prepare lookup table of sin and cos ( i*Lj )
  362.  * for required multiple angles
  363.  */
  364. int sscc( k, arg, n )
  365. int k;
  366. double arg;
  367. int n;
  368. {
  369. double cu, su, cv, sv, s;
  370. int i;
  371. double sin(), cos();
  372.  
  373. su = sin(arg);
  374. cu = cos(arg);
  375. ss[k][0] = su;            /* sin(L) */
  376. cc[k][0] = cu;            /* cos(L) */
  377. sv = 2.0*su*cu;
  378. cv = cu*cu - su*su;
  379. ss[k][1] = sv;            /* sin(2L) */
  380. cc[k][1] = cv;
  381.  
  382. for( i=2; i<n; i++ )
  383.     {
  384.     s =  su*cv + cu*sv;
  385.     cv = cu*cv - su*sv;
  386.     sv = s;
  387.     ss[k][i] = sv;        /* sin( i+1 L ) */
  388.     cc[k][i] = cv;
  389.     }
  390. return(0);
  391. }
  392.  
  393.